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Abstract We propose a finite volume scheme for convection-diffusion equations with nonlinear 
diffusion. Such equations arise in numerous physical contexts. We will particularly focus on the 
drift-diffusion system for semiconductors and the porous media equation. In these two cases, it is 
shown that the transient solution converges to a steady-state solution as t tends to infinity. 
The introduced scheme is an extension of the Scharfetter-Gummel scheme for nonlinear diffusion. It 
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scheme in the nondegenerate case. Finally, we present some numerical simulations applied to the 
two physical models introduced and we underline the efficiency of the scheme to preserve long-time 
behavior of the solutions. 
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1 Introduction 

In this article, our aim is to elaborate a finite volume scheme for convection-diffusion equations 
with nonlinear diffusion. The main objective of building such a scheme is to preserve steady-states 
in order to be able to apply it to physical models in which it has been proved that the solution 
converges to equilibrium in long time. In particular, this convergence can be observed in the drift- 
diffusion system for semiconductors as well as in the porous media equation. 

In this context, we will first present these two physical models - drift-diffusion system for semicon- 
ductors and porous media equation. Then, we will precise the general framework of our study in 
this article. 
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1.1 The drift-diffusion model for semiconductors 

The drift-diffusion system consists of two continuity equations for the electron density N(x, t) and 
the hole density P(x,t), as well as a Poisson equation for the electrostatic potential V(x,t), for 
t E K+ and x G R d . 

Let fi C K d [d > 1) be an open and bounded domain. The drift-diffusion system reads 

d t N - div(Vr(7V) - NVV) = on [2 x (0, T), 

d t P - div(Vr(P) + PVV) = on fix(0,T), (1) 
4V = iV-P-C on/2x(0,T), 

where C G L°°{f2) is the prescribed doping profile. 
The pressure has the form of a power law, 

r(s) — s 7 , 7 > 1. 

We supplement these equations with initial conditions N$(x) and Pq(x) and physically motivated 
boundary conditions: the boundary r = dfl is split into two parts P = r D U P and the boundary 
conditions are Dirichlet boundary conditions N, P and V on ohmic contacts r D and homogeneous 
Neumann boundary conditions on r(N), r(P) and V on insulating boundary segments r N . 
The large time behavior of the solutions to the nonlinear drift-diffusion model ([T]) has been studied 
by A. Jiingel in It is proved that the solution to the transient system converges to a solution 
of the thermal equilibrium state as t — > oo if the Dirichlet boundary conditions are in thermal equi- 
librium. The thermal equilibrium is a particular steady-state for which electron and hole currents, 
namely Vr(N) — NW and Vr(P) + PW , vanish. The existence of a thermal equilibrium has been 
studied in the case of a linear pressure by P. Markowich, C. Ringhofer and C. Schmeiser in [24jl23j . 
and in the nonlinear case by P. Markowich and A. Unterreiter in [25] . 
We introduce the enthalpy function h defined by 

h(s) = [ S ^-dT (2) 
h T 

and the generalized inverse g of h defined by 

„ M f^H*) if M0+)< S <co, 
9(S) -\ if a < h(0+). 

If the boundary conditions satisfy N, P > and 

h(N)-V = ot N and h(P) + V = ct P on r D , 

the thermal equilibrium is defined by 

N e "(x) =g(a N + V e «(xj) 1 P^{x) = g (a P - V^(xj) , x e Q, (3) 

while V eq satisfies the following elliptic problem 



f AV eq = g(a N + V eq ) -g{a P - V eq ) - C in Q, 
\V eq (x) = V(x) on r D , VV eq ■ n = on r N . 



(4) 
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The proof of the convergence to thermal equilibrium is based on an energy estimate with the control 
of the energy dissipation. More precisely, if we define 

H(s) = J h(r)dT, s > 0, (5) 

then we can introduce the deviation of the total energy (sum of the internal energies for the electron 
and hole densities and the energy due to the electrostatic potential) from the thermal equilibrium 
(see [20]) 

£(t) = J (h (N(t)) - H (N eq ) - h (N eq ) (N(t) - N eq ) + H (P{t)) - H (P eq ) 

-h (P eq ) (P(t) - P eq ) + \ |V (V(t) - V eq )\ 2 ^j dx, (6) 

and the energy dissipation 

At) = -J Q (N(t) Mh(N(t)) V(t))\ 2 + Pit) \V(HP{t)) + V(t))\ 2 ) dx. (7) 

Then the keypoint of the proof is the following estimate: 

< 8{t) + [ l{ T )dT < £(0). (8) 
Jo 



1.2 The porous media equation 

The flow of a gas in a d-dimensional porous medium is classically described by the Leibenzon-Muskat 
model, 

'd t v = Av~< onR d x(0,T), (] 
v(x,0) =v (x) onf d , [) 

where the function v represents the density of the gas in the porous medium and 7 > 1 is a physical 
constant. 

With a time-dependent scaling (see [7]), we transform ((9]) into the nonlinear Fokker-Planck equation 

d t u = div(ira + Vu 7 ) onR d x(0,T), , . 

u(a:,0)=uo(aO onR d . ^ ' 

It is proved in |7] that the unique stationary solution of pH|) is given by the Barenblatt-Pattle type 
formula 

/ ~_1 \ 1/(7-1) 

u*i(x) = (c 1 -l^-\x\ 2 ) , (11) 

where C\ is a constant such that u eq has the same mass as the initial data uq. 

Moreover, J. A. Carrillo and G. Toscani have proved in [7] the convergence of the solution u(x, t) of 

© to the Barenblatt-Pattle solution u eq (x) as t — > 00. As in the case of the drift-diffusion model, 
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the proof of the convergence to the Barenblatt-Pattle solution is based on an entropy estimate with 
the control of the entropy dissipation given by (|8|), where the relative entropy is defined by 



£{t) = j (H(u(t)) - H(u e «) + ^ (u(t) - dx, 

where H is defined by ([5]) and the entropy dissipation is given by 

d 



At) 



dt 



£[i) 



u(f) 



V h(u{t)) 



dx. 



(12) 



(13) 



1.3 Motivation 

Many numerical schemes have been proposed to approximate the solutions of nonlinear convection- 
diffusion equations. In particular, finite volume methods have been proved to be efficient in the 
case of degenerate parabolic equations (see [T51IT()] ). We also mention the combined finite volume- 
finite clement approach for nonlinear degenerate parabolic convection-diffusion-reaction equations 
analysed in [T7]. The definition of the so-called local Peclet upstream weighting numerical flux 
guarantees the stability of the scheme while reducing the excessive numerical diffusion added by 
the classical upwinding. 

On the other hand, there exists a wide literature on numerical schemes for the drift-diffusion equa- 
tions. It started with 1-D finite difference methods and the Scharfetter-Gummel scheme ([26]). In 
the linear pressure case (r(s) = s), a mixed exponential fitting finite element scheme has been 
successfully developed by F. Brezzi, L. Marini and P. Pietra in [SJH]. The adaptation of the mixed 
exponential fitting method to the nonlinear case has been developed by F. Arimburgo, C. Baiocchi, 
L. Marini in [5] and by A. Jiingel in [TJ5] for the one-dimensional problem, and by A. Jiingel and 
P. Pietra in [21] for the two-dimensional problem. Moreover, C. Chainais-Hillairet and Y.J. Peng 
proposed a finite volume scheme for the drift-diffusion equations in 1-D in [10] , which was extended 
in [9][TT| in the multidimensional case. C. Chainais-Hillairet and F. Filbet also introduced in [8] 
a finite- volume scheme preserving the large time behavior of the solutions of the nonlinear drift- 
diffusion model. 

Now to explain our approach, let us first recall some previous numerical results concerning the drift- 
diffusion system for semiconductors. The precise definitions of schemes considered will be presented 
in Section 2. We compare results obtained with three existing finite volume schemes: the classical 
upwind scheme proposed by C. Chainais-Hillairet and Y. J. Peng in [10], the Scharfetter-Gummel 
scheme introduced in |26j and the nonlinear upwind scheme studied in [8]. 

In Figure [TJ we present some results obtained in the case of a linear diffusion (r(s) = s). We rep- 
resent the relative energy £ and the dissipation of energy I obtained with the upwind flux and the 
Scharfetter-Gummel flux for a test case in one space dimension. We can observe a phenomenon of 
saturation of £ and T for the upwind flux. In addition, we clearly observe that the energy and its 
dissipation obtained with the Scharfetter-Gummel flux converge to zero when time goes to infinity, 
which means that densities N(t) and P(t) converge to the thermal equilibrium. It appears that the 
Scharfetter-Gummel flux is very efficient, but is only valid for linear diffusion. Moreover, we can 
emphasize that contrary to the upwind flux, the Scharfetter-Gummel flux preserves the thermal 
equilibrium. 

In Figure O we present numerical results obtained in the case of a nonlinear diffusion r(s) = s 2 . 
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Fig. 1 Linear case: relative energy £ n and dissipation I n for different schemes in log scale, with time step At = 10 
and space step Ax = 10~ 2 . 
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Fig. 2 Nonlinear case: relative energy £ n and dissipation I n for different schemes in log scale, with time step 
At = 5.f0 — 4 and space step Ax = f0~ 2 . 



We represent the relative energy £ and the dissipation I obtained with the classical upwind flux 
and with the nonlinear upwind flux for a test case in one dimension of space. We still observe a 
phenomenon of saturation of £ and X for the classical upwind flux. For the nonlinear flux, we clearly 
notice that the energy and its dissipation converge to zero when time goes to infinity. 
Looking at these results, it seems crucial that the numerical flux preserves the thermal equilibrium 
to obtain the consistency of the approximate solution in the long time asymptotic limit. 



Our aim is to propose a finite volume scheme for convection-diffusion equations with nonlinear 
diffusion. We will focus on preserving steady-states in order to obtain a satisfying long-time behavior 
of the approximate solution. The scheme proposed in [§] satisfies this property but because of the 
nonlinear discretization of the diffusive terms, it leads to solve a nonlinear system at each time 
step, even in the case of a linear diffusion. The idea is to extend the Scharfetter-Gummel scheme, 
which is only valid in the case of a linear diffusion, for convection-diffusion equations with nonlinear 
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diffusion, even in the degenerate case. Some extensions of this scheme have already been proposed. 
Indeed, R. Eymard, J. Fuhrmann and K. Gartner studied a scheme valid in the case where the 
convection and diffusion terms are nonlinear (see |13|). but their method leads to solve a nonlinear 
elliptic problem at each interface. A. Jiingel and P. Pietra proposed a scheme for the drift-diffusion 
model (see [T51I21) ). but it is not very satisfying to reflect the large-time behavior of the solutions. 



1.4 General framework 

We will now consider the following problem: 

d t u - div(Vr(u) - qu) = for {x, t) £ Q x (0, T), (14) 

with an initial condition 

u(x,0) = uq(x) for x £ O. (15) 

Moreover, we will consider Dirichlet-Neumann boundary conditions. The boundary dfl — T is split 
into two parts T = r D U JT and, if we denote by n the outward normal to r, the boundary 
conditions are Dirichlet boundary conditions on r D 

u(x, t) = u(x, t) for (a;, t) £ r D x (0, T), (16) 

and homogeneous Neumann boundary conditions on P : 

Vr(tt) ■n = 0onr A 'x(0, T). (17) 

Remark 1 We will construct the scheme and perform some numerical experiments in the case of 
Dirichlet-Neumann boundary conditions. However, for the analysis of the scheme, we will only 
consider the case of Dirichlet boundary conditions (dQ = r D = r). 

We suppose that the following hypotheses are fulfilled: 
(HI) fl is an open bounded connected subset of M. d , with d = 1,2 or 3, 

(H2) dfl = r D = r, u is the trace onfx (0, T) of a function, also denoted u, which is assumed 

to satisfy u £ H\f2 x (0, T)) fl L°°(f2 x (0,T)) and u > a.e., 
(H3) u £ L°°(f2) and u Q > a.e., 

(H4) r £ C 2 (R) is strictly increasing on ]0, +oo[, r(0) = r'(0) = 0, with r'(s) > c s 7_1 , 
(H5) q£ C 1 (Tl 1 R d ). 

H. Alt, S. Luckhaus and A. Visintin, as well as J. Carrillo, studied the existence and uniqueness of 
a weak solution to the problem ([T4|) - ([T7ll in [I] and [6] respectively. 

Definition 1 We say that u is a solution to the problem p4)) - ((T5|) - (jl6l) - p7|) if it verifies: 

u £ L°°(/2 x (0,T)), u-u £ L 2 (Q,T;H%(f2)) 
and for all ip £ V(Q x [0,T[), 

(ud t tp- V(r(u)) ■ VV' + uq- V^)da;dt+ / u(a:, 0) 0) da; = 0. (18) 
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The outline of the paper is the following. In Section 2, we construct the finite volume scheme. 
In Section 3, we prove the existence and uniqueness of the solution of the scheme and give some 
estimates on this solution. Then, thanks to these estimates, we prove in Section 4 the compactness 
of a family of approximate solutions. It yields the convergence (up to a subsequence) of the solution 
us of the scheme to a solution of (fT^)) - (fl"T)) when S goes to 0. In the last section, we present some 
numerical results that show the efficiency of the scheme. 



2 Presentation of the numerical scheme 



In this section, we present our new finite volume scheme for equation (|14[) and other existing 
schemes. We will then compare these schemes to our new one. 



2.1 Definition of the finite volume scheme 



We first define the space discretization of fl. A regular and admissible mesh of fi is given by a 
family T of control volumes (open and convex polygons in 2-D, polyhedra in 3-D), a family £ of 
edges in 2-D (faces in 3-D) and a family of points (xk)keT which satisfy Definition 5.1 in [15]. It 
implies that the straight line between two neighboring centers of cells (i/c, xl) is orthogonal to the 
edge a = K\L. 

In the set of edges £ , we distinguish the interior edges a £ £mt and the boundary edges a £ £ ext . 
Because of the Dirichlet-Neumann boundary conditions, we split £ ext into £ ex t — £^ xt U £^ xt where 
£® x t is the set of Dirichlet boundary edges and £^ xt is the set of Neumann boundary edges. For a 
control volume K £ T, we denote by £k the set of its edges, £mt,K the set of its interior edges, 
£® xt K the set of edges of K included in r D and £^ xt K the set of edges of K included in r N . 
The size of the mesh is defined by 

Ax = max(diam(X)). 

ifeT 

In the sequel, we denote by d the distance in R d and m the measure in R d or R d_1 . 
We note for all a £ £ 

_(d(x K ,x L ), for a £ £ inU a = K\L, 
a \d(x K ,a), for a £ £ ex t,K- 

For all er £ £, we define the transmissibility coefficient r CT ~ — - — -. 

For a £ £k, T^-K,a is the unit vector normal to a outward to K. 

We may now define the finite volume approximation of the equation (|14 |) -(|17 p . 

Let (T,£, (xk)kgt) be an admissible discretization of ft and let us define the time step At, Nt — 

E(T/At) and the increasing sequence (t n )o< n <N T , where t n = nAt, in order to get a space-time 

discretization T> of Q x (0,T). The size of the space-time discretization T> is defined by: 

S = max(Z\a;, At). 

First of all, the initial condition is discrctized by: 

U° K = -Ls [ u (x)dx, K£T. (19) 
m(-K) J K 
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In order to introduce the finite volume scheme, we also need to define the numerical boundary 
conditions: 

i r tn+1 r 

K l+1 = Atm{a) J n J u(s, t) ds dt, a e £g t , n>0. (20) 

We set 

Qk,<t = — rr / <l( x ) ■ ds{x), VK 6 T, Vo- 6 £ K . (21) 
ml 0- ) J a 

The finite volume scheme is obtained by integrating the equation (|14j) on each control volume and 
by using the divergence theorem. We choose a backward Euler discretization in time (in order to 
avoid a restriction on the time step of the form At = 0(Ax 2 )). Then the scheme on u is given by 
the following set of equations: 

jjn+X _ jjn 

m{K) At + £ ^ 1 = < ( 22 ) 

a££ K 

where the numerical flux J~^ + ^ is an approximation of — / (Vr(u) — qu) ■ n#- jCr which remains to 

J a 

be defined. 



2.2 Definition of the numerical flux 
2.2.1 Existing schemes 

We presented in introduction some numerical results obtained with different choices of numerical 
fluxes for the drift-diffusion system. We are now going to define precisely these fluxes. 

The classical upwind flux. This flux was studied in [15] for a scalar convection-diffusion equation. 
It is valid both in the case of a linear diffusion and in the case of a nonlinear diffusion. The diffusion 
term is discrctized classically by using a two-points flux and the convection term is discretized with 
the upwind flux, whose origin can be traced back to the work of R. Courant, E. Isaacson and M. 
Rees [12J . This flux was then used for the drift-diffusion system for semiconductors in [10] and [9] 
[TT] in 1-D and in 2-D respectively. The definition of this flux is 

r« (r (t/r 1 ) - r (U2 +1 ) + d a U a Ul +1 qj^U? 1 )) , ^ = K\L G £ int , K , 
F5& = {rjr - r (U^) + d a U a U^ - q^A) - ^ 6 (23) 



Va e £ N 



ext.Ki 



where s + = max(s, 0) and s~ = max(— s, 0) are the positive and negative parts of a real number s. 
The upwind flux with nonlinear discretization of the diffusion term. This flux was intro- 
duced in [S] in the context of the drift-diffusion system for semiconductors. The idea is to write 

the flux — / (Vr(u) — qu) ■ n^ CT as — / (uVh(u) — qu) ■ n.K,a, where h is the enthalpy function 

J a J a 
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defined by ([2]). The flux is then defined with a standard upwinding for the convective term and a 
nonlinear approximation for the diffusive term: 

-pn+l _ 

-r a (min (U%+\ U n L +1 ) Dh (U^) K a + d a U, a U n K +1 - q^U^X) , Va = K\L, 

-r a (min U^) Dh (t/" +1 )^ + d a (<?+ ff £/£ +1 - q^US +1 )) , Va E £g 

0, Va € £, w 

where for a given function /, Df(U)K.a is defined by 



ext,Ki 
■N 

ext,K' 



f(U L ) - f(U K ), if a = K\L e £ 



K.int i 



Df{U) K . a = { f(U a ) - f(U K ), tfae£ D 
0, if a € £ 



K,ext ' 
N 

K \ext ' 



This flux preserves the thermal equilibrium and it is proved that the numerical solution converges 
to this equilibrium when time goes to infinity. 

The Scharfetter-Gummel flux. This flux is widely used in the semiconductors framework in 
the case of a linear diffusion, namely r(s) = s. It has been proposed by D.L. Scharfetter and 
H.K. Gummel in [2 6) for the numerical approximation of the one-dimensional drift-diffusion model. 
We also refer to the work of A.M. Il'in [18] , where the same kind of flux was introduced for one- 
dimensional finite-difference schemes. The Scharfetter-Gummel flux preserves steady-state, and is 
second order accurate in space (see [22 ). It is defined by: 

T a {B{-d a q K . a )Ul +1 - B(d a q K . a )U2 +1 ) , Va = K\L e £ KMU 
F£* = <! T - (B(-d„q K ^U n K +1 - B(d a q K ,M +x ) , Va e £° 



o, Va e £ N 



K,exf 



where B is the Bernoulli function defined by 



BW = — fon/0, B(0) = 1. 
e x — 1 



2.2.2 Extension of the Scharfetter-Gummel flux 



Now we will extend the Scharfetter-Gummel flux to the case of a nonlinear diffusion. Firstly, if we 
consider the linear case with a viscosity coefficient e > 0, namely 

d t u - div(eVu - qu) = for (x, t) e H x (0, T), 

then the Scharfetter-Gummel flux is defined by: 

r^^r^^B^^^U^-B^f^^U^^ W = K\Le£ mt ., K . (24) 

Using the following properties of the Bernoulli function: 

B(s) — > and B(s) s, 

s^ + oo — oo 
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it is clear that if e tends to zero, this flux degenerates into the classical upwind flux for the transport 
equation dtu — div(qu) = 0: 



^K.cr 



Tjn+l 
u K 



n+1 
L 



V(7 = K\L G £ 



int.K ■ 



(25) 



Now considering a nonlinear diffusion, we can write Vr(u) as r'(a)Vii. We denote by drK,a an 
approximation of r'(u) at the interface a G £k, which will be defined later. We consider this term 
as a viscosity coefficient and then, using (|24p . we extend the Scharfetter-Gummel flux by defining: 



x-n+l 
^K,a 



T a dr K ,< 

T a dr K> , 
0. 



-d a qK,cr 
dr K ,G 

-d a qK,cr 
dr 



U K +1 - B 



U K +1 - B 



d a qK,a 

dr K ,a- 

d a qK,a 

dr 



ul+ i 



u, 



n+1 



, Ver = K\L G £int,K, 



(26) 



In the degenerate case, drK,a can vanish and then this flux degenerates into the upwind flux (|25jl. 
Now it remains to define drK a ■ 



Definition of drx,a- A first possibility is to take the value of r' at the average of Uk and U a : 

'Uk + Ul 



dr K ,a = 



u K + u a 



Vcr = K\L G £, 



int.K i 



g e° 



(27) 



This choice is quite close to the one of A. Jiingel and P. Pietra (see [191121] ). However, considering 
the numerical results presented in the introduction, it seems important that the numerical flux 
preserves the equilibrium. Therefore, we define the function dr as follows: for a, b G R+, 



dr(a, b) 



h{b) - h{a) 

log(6) - log(a) 
'a + b 



and we set for all K G T 



A', i 



if afc > and a ^ b, 
elsewhere, 



(dr(U K ,U L ), ioi a = K\L E £ K ,int, 
\dr(U K ,U a ), forae£° ext . 



(28) 



(29) 



Remark 2 Let K <E T and cr G £a". We assume that drK, a is defined by (1291) in (|26|) and that 

Uk > and £7 CT > 0. If d a q K ,a = Dh(U) K ,a, then = 0. 

Indeed, 



^a,o- = r a dr K .a \ B - 



Dh(U) K , 



dr 



Uk-B 



(Dh{U) K , 



T a Dh{U) K ,, 



V dr A, 

ex P — 3 (7a - \ 

V « r A,<x / 



0* 



exp 



r ph(u) K ,a 
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But using the definition (|28[1 of dr, we obtain 

/ Dh(U) K ,a \ = 

and then Tk,o = 0. Thus the scheme preserves this type of steady-state. 
Time discretization. We choose an explicit expression of drx,a'- 

_(dr(U%,UZ), for a = K \L e £ K ,int, ( ^ 

Thus we obtain a scheme which leads only to solve a linear system of equations at each time step. 
To sum up, our extension of the Scharfetter-Gummel flux is defined by 



•r-n+l 



~d a qK,cr 


u n+1 


— B 


d G qK,<j 


dr n 


dr n 
\ K,a 


—d a qK. a 


K 


- B 


d a qK,(j 


dr n 

K,a 1 


dr n 



(31) 



where cir^- is defined by (j30l) . This flux preserves the equilibrium. 
2.3 Consistency of the numerical flux 

Lemma 1 Let a, tel, a,b > 0. TTien there exists rj € [min(a, 6), max(a, 6)] smc/i £/ia£ 

dr(a, 6) = r'(rf). 

Proof The result is clear if ab — or a = b. Let us suppose that ab > and a < 6 (the proof is the 
same if a > b). If we consider the change of variables x = log(o) and y = log(6), we obtain 

= h(exp(y)) - fo(exp(a;)) 
y-x 

and using Taylor's formula, there exists 9 G [x, y] such that 

dr(a,b) = exp(9)h' (exp(9)) = r'(exp(9)) (using the definition of h). 
Finally, there exists r\ = exp(9) S [a, b] such that 

dr(a, b) = r (rj). 
Remark 3 The flux ([31]) can also be written as 

U K +1 + m(o-)q K ,a „„ +u ( d*q_K. 

2dr 



rf) = m[o)q K ,„ K Z ° t^J^L coth ^ {U n + l un + iy (32) 



K.i 



The first term is a centred discretization of the convective part. The second term is consistent with 
the diffusive part of equation (|T4|) . since coth(x) ~ — . 
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3 Properties of the scheme 

3.1 Well-posedness of the scheme 

The following proposition gives the existence and uniqueness result of the solution to the scheme 
defined by ([l9] ) -([2D ll -(|22 |) - ([3T|l and an L°°-estimate on this solution. 

Proposition 1 Let us assume hypotheses (H1)-(H5). Let T> be an admissible discretization of Q x 
(0,T). Then there exists a unique solution {U^,K G T, < n < Nt} to the scheme i!9\) - [2l^) - 
fUj-fHP, with [/£ > for all K G T and <n < N T . 
Moreover, if we suppose that the two following assumptions are fulfilled: 

(H6) div(q) = 0, 

(HI) there exist two constants m > and M > such that m <u,Uq < M, 
then we have 

0<m<U£< M, VK 6 T, Vn > 0. (33) 



Proof At each time step, the scheme (|T9l) - (|20p - (p2)) -([3~T |) leads to a system of card(T) linear equations 
on U n+1 — (Uk +1 )keT which can be written: 

A n U n+1 = S n , 



where : 

• A" is the matrix defined by 



A n K , L = - r c dr n K ^B I VLeT such that a = K\L 6 S int , K ; 



dr n 



At 



• S n = ( ^t^Uk I + Tb n , with 



KeT 



if K G T such that m(dK n f) = 0, 



Tb\ = 



V r CT drJ CT S ( ) C/; l+1 if K G T such that m(8K fl T) / 0. 



The diagonal terms of A n are positive and the offdiagonal terms are nonnegative (since B(x) > 
for all x G M. and dr\ > for all KeT, for all er G £k)- Moreover, since dr\ — dr*l a and 
qK,a = —q_L,a for all a = K\L G £int, we have for all LeT: 



KeT 



and then A" is strictly diagonally dominant with respect to the columns. A n is then an M-matrix 
so A n is invertible, which gives existence and uniqueness of the solution of the scheme. Moreover, 
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{A" 1 )- 1 > and since U% > for all K G T (using (H3)) and i7™ +1 > for all n > 0, for all 

a G S^xt (using (H2)), it is easy to prove by induction that 11% > for all K 6 T, for all n > 0. 

Now, we suppose that (H6) and (H7) are fulfilled. We prove that U% < M for all K G T, for all 

n > by induction. Thanks to hypothesis (H7), we have clearly Ug- < M for all K e T. 

Let us suppose that U% < M VK G T. We want to prove U^ +1 < M \/K G T. 

Let us define M = (M, M) T G R card ( r ). Since A n is an M- matrix, we have (A 11 )- 1 > and then 

it suffices to prove that A n (U n+1 - M) < 0. 

We first compute -A"M. Using the following property of the Bernoulli function: 

B{x) - B(-x) = -x Vx G M, (34) 

we obtain that for all K G T, 



— -t- 2^ mi ";^^ -t- 2^ T ° ar k,<r D i - 

Then we compute A n (U n+1 - M): for all K G T 



-M m(tr) 9jr , a - Af £ r a ir\ a B I 

cr££ int ,K <T&ER. 4 u- \ 



dr n 

K,(7 



By induction hypothesis, the first term is nonpositive. Moreover, using hypothesis (H7) and the 
property (|34|). we obtain 

(A n (U n+1 -M)) K < -M H°)<lK, a -M J2 M?)Qk,c 

< -M m(o)q K ,*. 

<t&£k 

However, using hypothesis (H6) and the definition of qx,a we get 

ra(a)q K , a = ^ / 1 ' n ^ rfs = / div (?) = °> 
and then (A" (U n+1 - M)) K < for all K <eT- 

So we have A™ (C/" +1 - M) < 0, therefore we deduce that U n+1 - M < 0, hence Up" 1 < M \/K 
and we can show by the same way that Up" 1 > m VK. 

Remark 4 In the case of the drift-diffusion system for semiconductors, the hypothesis (H6) is not 
fulfilled (AV ^ 0). Nevertheless, if we assume that 

— the doping profile C is equal to 0, 

— there exist two constants m > and M > such that m < N , Nq, P, Pq < M, 

— MAt < 1, 
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then we have, using the same kind of proof as in [5], 

< to < N£ < M, VK G T, Vn > 0, 
< m < P£ < M, G T, Vn > 0. 

Definition 2 Let T> be an admissible discretization of fl x (0,T). The approximate solution to 
the problem (fTl)) - (fTI^ - ([TB| - (fTTl) associated to the discretization V is defined as piecewise constant 
function by: 

u s (x,t) = Ul +l , V(x,t) eK x[t n ,t n+1 [, (35) 
where {U%, K 6 T, < n < N T } is the unique solution to the scheme (fI5 ]) -([2l? ]) -P ^ -(|3T ]) , 



3.2 Discrete L 2 (OjT;!! 1 ) estimate on us 

In this section, we prove a discrete L 2 (0, T; H 1 ) estimate on us in the nondegenerate case, which 
leads to compactness and convergence results. 

For a piecewise constant function vs defined by vg(x,t) = v 7 ^ 1 for (x,t) G K x and 
V 5 (7,f) = < +1 for (7,t) G cr x [*V n+1 [> we define 



A'r 



/ 



TO 



ll,2> = E^ 



n=0 



\a=K\L 



\ 



J 



Proposition 2 Let assume (H1)-(H7) are satisfied. Let us be defined by the scheme fl9\) -$20 \) - 
(EM>-(3^ and fjf^. 

There exists D\ > only depending on r, q, Uq, u, fl and T such that 



\us\\\.v < Di. 



(36) 



Proof We follow the proof of Lemma 4.2 in [T3]. Throughout this proof, Di denotes constants which 
depend only on r, q, uq, u, Q and T. We set 



T7"+l 



and 



We multiply the scheme 



w 



n+1 - UI+ 1 - U n x\ VK G T, Vn G N. 



A" 



by Atw 1 ^ 1 and we sum over n and K. We obtain A + B = 0, where: 



AT-r 



A = E E m (*) - ^ 

n=0KeT 

AT T 

B = E^E E ^r>* +1 



a: J w at i 



n=0 KeTaeSi, 
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Estimate of A. This term is treated in [T3]. We get: 

A> -^\\uo-u(.,0)\\h W -nd t u\\ LHnx{ o,T))\M-m\ = -D 2 . (37) 

Estimate of B. A discrete integration by parts yields (using that u>™ +1 = for all a G £® xt 
and for all n > 0): 



B = T. At E -^K +1 -< +1 ) + E^E E t k) R+ 1 - > 

(7—K\L 

which delivers B = B' — B, with: 

B ' = j2*t E ^ - u l +1 ) + E ^ E E ^ ( u k +1 u « +1 ) , 

n=0 cre£i„t "=0 KeTa££ D . „ 

(7 — iv |L 

n=0 <T6£ ml n=0 KeTae£ D , w 

a—K\L 

Estimate of B. Using the expression ([32]) of J 7 ^^ , we have B — B\ -\- B2 with 

5i = E ^ E + u n - uT 1 ) 

n — &££int 

a=K\L 

5.=e^e -th (W) - lt-h) 1 - c/r x ) 

cr—K\L 

+X>E E ^^coth^)^ 1 -^ 1 )^ 1 -^ 



The term £?! is treated like in [13], which leads to 

|S 1 |<M[|q[[ 00 [[u5||i^dm(^)=- D 3 
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We apply Young's inequality for B^'. for any a > 0, we have 

+?X>£ £ -(*U 2 (|fW|^))Vr-^') 2 

n=0 K£To-E£V t K \ K ' a \ K ' a / / 

+ 2^INII?,z>- 

By the hypothesis (H4), we have inf r'(s) > 0. Then, using Lemma HJ the L°° estimate on us 

s£[m,M] 

([33]) and the hypothesis (H5), we have 

dgqK,a . ||q||oodiam(l?) ^ M <r w 

Moreover, since a; i— >• x coth(a;) is continuous on K, we obtain 

(^«*> (^)) 2 < v„ e N, V* E T, V. e &. 

Thus we can bound B: 

\B\ < D 3 + ^£> 4 | sup r'(s)| ||uj||? 1J + (38) 



Sg[TO,M] 



2a 1 



Estimate of B' . First, using the expression (|3^|) of the flux and LemmaQ] we have for all n > 0, 
for all K £ T and for all <r = £ £int,K 



(r} Kttr ) cotn \ [U K U L ) 



Then, since x coth(x) > 1 for all x £ K, we get: 

m(a)q K ,a 

•uc. ! ■ - i ' ' / i i -.- v. T 

sS[m,M] 

We obtain the same type of inequality for J"^ 1 (U^ +1 - C/™ +1 ). Thus we get 

N< 



- C/r 1 ) > ((^ +1 ) 2 - (U2 +1 ) 2 ) + T„ s inf r'(s) (U^ - U^f 



^^/MiMi^+E^ E ^^((^ +1 ) 2 -K +1 ) 2 ) 

n=0 cre£i„ t 

f>E E ^^((^ +1 ) 2 -(^ +1 ) 2 ). 
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Through integrating by parts and using the hypothesis (H6), we get 

m{a)q Ki<T ( (Trn+u 2 , rm+ n2 



a=K\L 

m(a)q K , 



L 



5>E E ^f^((^ +1 ) 2 -(^ +1 ) 2 ) 
X> E E \ I ■ n ^ ds ^ (^ n+1 ) 2 = 

n=0 KGT^afO 



and then 



B'> inf r'OOIKH^-A;. (39) 
Conclusion. Using ^4 + £> = and estimates (|3T|) . ([35]) and (|39p . we finally get for any a > 0: 

W /'(s)-^4 sup r'(s)) )\\us\\l v <D 2 + D 3 + D 5 + ^-\\u s \\ 2 lv , 
se[m,M] 2 ys£[ m ,M] J I 2a 

2 inf r'(s) 

s£[m,iW] II 9 

thus for a < we obtain \\us\\l j, < D%. 

D 4 sup r'(s) 

\s<E[m,M] I 



4 Convergence 

In this section, we prove the convergence of the approximate solution us to a weak solution u 
of the problem (fT4l) - (pT5|) - (P~^|l - (|T7 |) . Our first goal is to prove the strong compactness of (us)s>q 
in L 2 (J?x]0,T[). It comes from the criterion of strong compactness of a sequence by using esti- 
mates (|33| and ((36|) . Then, we will prove the weak compactness in L 2 (f2x]0,T[) of an approximate 
gradient. Finally, we will show the convergence of the scheme. 



4.1 Compactness of the approximate solution 

The following Lemma is a classical consequence of Proposition [5] and estimates of time translation 
for us obtained from the scheme (fT9^) - ([2H| - P2^) - ([3T|) . The proof is similar to those of Lemma 4.3 
and Lemma 4.7 in [15] . 

Lemma 2 (Space and time translate estimates) We suppose (H1)-(H7). Let V be an admis- 
sible discretization of fl x (0,T). Let us be defined by the scheme H19\) - l£7jj) - i22\) - Wf]) and by \35\ ). 
Let u be defined by its — us a.e. on Q x (0,T) and us = a.e. on M. d+1 \ fl x (0,T). 
Then we get the existence of M 2 > 0, only depending on Q, T, r, q, it , u and not on T> such that 

I [ (us(x + n,t) -u s (x,t)) 2 dxdt < M a |»?|(|»j| +4<5), Vi; 6 f d , (40) 
Jo Jn 
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T p 

/ (u s (x,t + t) -u 5 (x,t)f dxdt < M 2 \t\, Vrel. (41) 
o J n 

Now, we define an approximation V s u$ of the gradient of u. Therefore, we will define a dual 
mesh. For K € T and a G f^-, we define Tr- iIT as follows: 

— if a = K\L G £mt,K, then Tr-.o- is the cell whose vertices are xk, xl and those of a — K\L, 

— if a 6 £ ex t,K, then Tx,o- is the cell whose vertices are a;^ and those of <r. 

See |llj for an example of construction of Tk,<j- Then ({Tk.o) defines a partition of Q. 
The approximation \7 s u$ is a piecewise function defined in Q x (0, T) by: 



V a u 5 (a:,t) 



' n(rT| (C/r +1 - "k, if (M) 6 ?K,ff x [fV" +1 [, a = X|L, 



-^^-(UZ +1 -Ul +1 )n K , a i{(x ) t)eT Ktt7 x[t n ,t n+1 [, ae£ ext , K . 



Proposition 3 We suppose (H1)-(H7). 

There exist subsequences of (ug)s>o and (V us)s>o, still denoted (us)s>o and (V us)s>o, and a 
function u 6 L°°(0,T; H l (n)) such that 

u$ — y u in L 2 (f2x]0,T[) strongly, as S — > 0, 
V 5 u s -± V« in {L 2 (f2x]0,T[)) d weakly, as S -> 0. 



Proof Using estimates (|4CH - (j4T|) and applying the Riesz-Frechet-Kolmogorov criterion of strong 
compactness [5], we obtain the first part of this Proposition. The result concerning V us is proved 
in [9]. 



4.2 Convergence of the scheme 

Now it remains to prove that the function u defined in Proposition^ satisfies Definition [T] of a weak 
solution. The main difficulty in proving this comes from the fact that the diffusive and convective 
terms are put together in the Scharfetter-Gummcl flux. 

Theorem 1 Assume (H1)-(H7) hold. Then the function u defined in Proposition^ satisfies the 
equation \lJ$ -ll5 }) - hlb}) - \17^ in the sense of H18\) and the boundary condition u—u G _L°°(0, T; H^fi)). 

Proof Let G V(tt x [0,T[) be a test function and i/>£ = tp[x K ,t n ) for all K e T and n > 0. We 
suppose that S > is small enough such that Supp(-0) C {x E (1; d(x,T) > 6} x [0, (Nt — 1) zit [. 
Let us define an approximate gradient of ip by 

m ( cr ) r.i.n „/.n \ „ -,c /_ x\ ^ „ i>n j-n+lr 



S.h(* t\ - J m ( 2 >f.«r. 



v>om) = 



m(er) 
m(T^cr' 



(VE - njf^ if (as, t) G Tk, ct x [i n , T +1 [, a = K\L, 

(C " 1%) n K,a if (I, t) G Tftfj X [<", T+ 1 [, a G 
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We get from [T3] that {V s i())s>o weakly converges to Vi/> in {L 2 {Q x (0,T))) d as 8 goes to zero. 
Let us introduce the following notations: 

Bio(S) =— I / / u$(x,t) dttp(x,t) dx dt + / us(x, 0) tp(x, 0) dx , 



Wo Jn Jn I 

B 20 (S) = [ [ r'{u s {x,t - At))\7 s u s (x,t) ■ Vif)(x,t) dx dt, 
Jo Jn 

B 30 (S) = - f [ u s (x,t)q(x) ■ \J s ip(x,t)dxdt, 
Jo Jn 

and 

s(S) = -B w (5) - B 20 (6) - B W (S). 
Multiplying the scheme by Atip^ and summing through K and n, we obtain 

B^S) + B 2 (S) + B 3 {S) = 0, 



where 



= E E < K ) ( u k +1 - u$) 

n=0KeT 

B 2 (S) = - £ At E E ^^T^ -th ( |^ ) (E^i - ^) *S 

n=0 KeTae£ K \ Zar K.c 



K • 



Nt TT n+1 4- 77™+! 

n=0 K£T<j££k 

From the strong convergence of the sequence (us)s>o to u in L 2 (i?x]0, T[), it is clear using the time 
translate estimate (|4ip that there exists a subsequence of (ug)s>o, still denoted by (us)s>o, such 
that 

us( ■ , ■ - At) — >u in L 2 {Hx}0, T[) strongly as S -> 0, 

where it 6 L ca {Q,T;H 1 (Q)) is defined in Proposition [3j Moreover, thanks to hypothesis (H4), we 
have r' G C 1 (M), and using the L°°-estimate ([55)) we obtain that 

r'(u s ( • , • - /It)) — ► r'(u) in L 2 (/2x]0, T[) strongly as <5 -> 0. 

Finally using this strong convergence and the weak convergence of the sequences (\7 s us)s>o to Vw 
and (VV)<5>o to Vip in (L 2 (/2x]0, T[)) d , it is easy to see that 

e(S) — > [ [ (u(x, i) dfip ~ r'(u(x, t)) Vu(x, i) ■ Vip + u(x, t) q(x) ■ Vifj) dx dt 
Jo Jn 

u(x, 0) ip(x, 0) dx, as 5 —> 0. 

Therefore, it suffices to prove that e(S) — > as 5 — >• and to this end we are going to prove that 
e(S) + B^S) + B 2 {5) + B 3 (S) — > as S -> 0. 
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Estimate of B\{8) — B\q{8). This term is discussed for example in [5] (Theorem 5.2) and it is 
proved that: 

\BxiS) - B 10 (S)\ < [(T + l)m(tf)M||V>||c»(nx(o,T))] S — ► as 5 -> 0. 
Estimate of 62(8) — B2o(S). Using a discrete integration by parts, we write 

b 2 (s) = £ At -th K +1 u^) m - r K ). 

a=K\L 

Then we rewrite B 2 (S) = B 2 i(S) + B 22 (8) + B 23 (<5), with 

B 21 {8) = Y / A t W(u%)(u2 +i -ui +i ){r L -r K ), 

o=K\L 



N T / 

B 2 2(6) = Y, At E t A 

n=0 a££ int \ 



2dr%„ \ 2dr K 



cr=K\L 

B 23 (5) = Y, At E r*(dr$,a-r'(U$))(Up 1 -Up 1 )y>l-n)- 

o=K\L 

Using the definition of us and V s us, we rewrite B2o{5) as i?2io(<5) + -B22o(<5) with: 
B 210 (S) = Y1 E r'm)^K(U2 +1 -U^ +1 ) / / WOr^.n^dzdi, 

n=0 CTG £„ lt m U^J ./*> Mr.- 

<t=K\L 



Now we prove that B 2 i(S) - B 2 w{5) -4 as <5 -4 and B 2 2{8), B 2 3(8), B 2 2o(8) -4 as <5 -4 0. 
Estimate of B 21 (<5) — i^io^)- We have 



Bai(*) - Saio(«) = E E m W r '(^) 

n—0 <j££int 



Since the straight line xk%l is orthogonal to the edge K\L, we have — xk = dcrn_ff !(T and 
then from the regularity of ifi, 



W^(x K ,t n )-n K>(T + 0(Ax) 



= Vtl>(x,t) ■ n K . a + 0(8), V(x,t) G T^ )CT x . 
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Then by taking the mean value over Tr- j(T , there exists Dq > depending only on tp such that 



n-n 



m(Tk l0 -) Jtk, 



VV> • nj<- a dx ) dt 



and then 



\B 21 (6) - B 210 (6)\ <6D 6 sup r'(s)J2^t £ m(cr) | C/£ +1 - C/£ +1 1 . 

s£[m,M] n=Q CTg£irit 

Since the straight line xk%l is orthogonal to the edge a — K\L for all er g £i n t,K and the mesh 
is regular, there is a constant Dj > depending only on the dimension of the domain and the 
geometry of T such that m(a)d a < Dtm(TK,a) for all K G T , all a E £ e xt.K and then using the 
Cauchy-Schwarz inequality and the L 2 (0,T; H 1 ) estimate we obtain 



\B 21 (8) - B 2W (8)\ < 8D e sup r'{s)y / D 1 TD 7 m(Q) — > as 8 ->■ 0. 

s£[m,M] 

Estimate of B 22 {8). Since x H> xcoth(a;) is a 1-Lipschitz continuous function and is equal to 
1 in 0, we have 



N T 



m(cr) 



qKAK +1 -u n K +i \\rL-r K \ 



\B 22 (8)\<^At , 

71=0 cre£ int 

< 25||q|U J2^t T ° K +1 ~ U^ +1 \ Wl - $k\ , since d a < 28. 



n—Q o^&int 

Then using the Cauchy-Schwarz inequality, the regularity of ip and the £ 2 (0, T; H 1 ) estimate (|36l) . 
there exists Ds > only depending on T and i? such that: 

\B 22 (5)\ < ^Hqlloo^sll^llci \/D~i — > as 5 -> 0. 

Estimate of -623(1$) • Using Lemma [T] and hypothesis (H4), we have 

\dr n K , a - r'(t/£)| < sup |r»| \U£ - , W e £ int) a = tf|L. 

s£[m,A/] 

Using the regularity of "0 and the Cauchy-Schwarz inequality, we obtain 

\B 23 (5)\<5 sup \r"(s)\Mc^At £ r CT \U£ - U%\ |C7£ +1 - Z7£ +1 | , 

and then using the L 2 (0,T; H 1 ) estimate (l36l) . we get 

|-B 23 (<5)| < * sup |r"(s)|||^[| C i£)i — > as 5 ->■ 0. 

se[m,M] 

Estimate of -B220 (<!))• We obtain the same type of estimate as for B 2 ^{8): 



\B 220 {8)\ < 28 sup |r"(s)|||^|| c iDi — > as 8 -> 0. 

se[m,M] 
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Estimate of S3 (6) — B 30 (S). Using a discrete integration by parts, we obtain 



Nt _ JT n+1 -I- TT n+1 

and then we rewrite B 3 (5) as -831(5) + B 32 (S), with 

N T JJ n+1 - JJ n+1 

B 31 (S) = -^At J2 m(a)q K , a — 2 K - Vk) 

n=0 crESint 
N T 

B 32 (s) = -j^At h*)<ik^uz +1 (r L - r K ) ■ 

n=0 <rE£int 



Using the definition of V 5 V; we get 

Nt 



B 30 (S) 



T .- r r 



Wl ~ ^k) <l(x) ■ n K,a dx dt, 



n=0ae£, 

which gives, using the definition of us, B 30 (S) = B 3 i (S) + B 320 (S), where 

B 3W (S) = -X> E m ( ff ) ( U l +1 U k +1 ) W£ - ^k) -7^~, I q(«) • ik, dx, 

n=0 <r€£ 4 „ t ^-LK,*) JT K ,„nL 

B 320 (S) = - E E m (^)^ +1 - V£) -77^ / q(*) • <fc. 



Now we prove that B 32 (5) - B 320 (S) -> as <5 -> and S 3 i(5), S 3 i (5) -> as <5 -> 0. 
Using the regularity of q, there exists Dg > which does not depend on S such that 



■7-r / q(z) • nff i(I ds(x) -j; — - / q(x) • n K , CT 

Then we can estimate B 32 (S) — B 320 (S): 



< D 9 S. 



N T 

\B 32 (S)~B 32a (5)\<SD 9 AlY,At m(a)|VE-^l 

< 8D s D 9 M\\iP\\ c i ^/D 7 ui{Q) — > as S -> 0. 



Moreover, we have 



Nt 



|B 3 i(5)|<5||q||ocE^ E ^I^-^IIV'E-V'S'I 

n—0 &E£in.t 

< 5\\q\\ 00 \\iP\\ el D 8 ,/D~ 1 — > as <5 -> 0. 
We obtain in the same way that i?3io(<5) — > as (5 — ► 0. 
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Hence u satisfies 




(u(x, t) d t ip(x, t) + r'(u(x, t)) Vu(x, t) ■ V<ip(x, t) + u(x, t) q(x) ■ Vip(x, t)) dx dt 




and then 




Jo Jn 




J a 



It remains to show that u - u G L°°(0, T; Hq([2)). This proof is based on the L 2 (0,T; H 1 ^)) 
estimate (j3"u]) and is similar to the one of Theorem 5.1 in [3]. 

5 Numerical simulations 

5.1 Order of convergence 

We consider the following one dimensional test case, picked in the paper of R. Eymard, J. Fuhrmann 
and K. Gartner [13]. We look at the case where, in flU]) we take fl = (0, 1), T = 0.004, r : s ^ s 2 , 
q = 100, in (fl"5)) we take uq = and in P^)l we take, for v — 200, 



The time step is taken equal to At — 10 -8 to study the order of convergence with respect to the 
spatial step size Ax. In Tables [1] and [2] we compare the order of convergence in L°° and L 2 norms of 
the scheme (fll?|) - (f2T)|) - (f2"2"l) defined on one hand with the classical upwind flux (|2"5|) and on the other 
hand with the Scharfetter-Gummel extended flux (|31[) . We obtain the same order of convergence 
as in [13] . Moreover, it appears that even if we are in a degenerate case, the Scharfetter-Gummel 
extended scheme is more accurate than the classical upwind scheme. 

5.2 Large time behavior 

5.2.1 The drift- diffusion system for semiconductors 



u(0,t) = (v- q)vt/2 




fori<l/u, 
(v — q)(vt — l)/2 otherwise. 



The unique weak solution of this problem is then given by 




We may define the finite volume approximation of the drift-diffusion system ([T]) . Initial and bound- 
ary conditions are approximated by (|19p and (|20p. The doping profile is approximated by (Ck)k^t 
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Ax(j) 


||« - Wj||ioo 


Order 


II" - Wj||z,oo 


Order 






Upwind 




SG extended 







2.5. 


1.110 




2.137.10- 1 




1 


1.25.10- 2 


7.237.10- 1 


0.62 


1.107.10- 1 


0.95 


2 


6.3. 10" 3 


4.485. 10" 1 


0.69 


5.631.10" 2 


0.98 


3 


3.1. 1CT 3 


2.685. 10" 1 


0.74 


2.84.10- 2 


0.99 


4 


1.6. 10~ 3 


1.568. lO" 1 


0.78 


1.426.10- 2 


1 


5 


8.10" 4 


9.10~ 2 


0.80 


7.15. 10~ 3 


1 



Table 1 Experimental order of convergence in L°° norm for spatial step sizes Ax(j) = of the classical upwind 

scheme and of the Scharfettcr-Gummel extended scheme. 
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Ax(j) 


\\U - U 6 \\ L 2 

Upwind 


Order 


\\U - U S \\ L 2 

SG extended 


Order 





2.5. 10~ 2 


3.336. 10- 1 




4.806.10-2 




1 


1.25.10-2 


1.852.10" 1 


0.85 


1.642.10-2 


1.55 


2 


6.3. 10" 3 


9.911.10" 2 


0.9 


5.695.10" 3 


1.53 


3 


3.1. 10~ 3 


5.182.10" 2 


0.94 


2.10" 3 


1.51 


4 


1.6. 10~ 3 


2.669.10" 2 


0.96 


7.142. 10" 4 


1.49 


5 


8.10" 4 


1.361.10-2 


0.97 


2.695.10" 4 


1.41 



Table 2 Experimental order of convergence in L 2 norm for spatial step sizes Ax(j) = ^ J+2 °f * ne classical upwind 
scheme and of the Scharfettcr-Gummel extended scheme. 



by taking the mean value of C on each volume K. The scheme for the system ([1} is given by: 

jyn+l _ iyra 

™(K) K At K + E T K + a = 0, V/v 6 T, Vn > 0, 

ae£K 

( pn+1 _ pn 

m(K) K — K + £ Vk) = °. V ^ e r ' V " ^ °> 

_ ae£K 

. J2*es K ^DV^ = m(K) (JV£ - Pg - C K ) , G T,Vn > 0, 

where 

/ / —DV n \ f DV n \ \ 

and 

/ / DV n \ ( —DV n \ \ 

We compute an approximation (N 1 ^ , P^ 1 ,V^)kgT of the thermal equilibrium (N eq , P eq ,V eq ) de- 
fined by ©-(jlj with the finite volume scheme proposed by C. Chainais-Hillairet and F. Filbet in 

Then we introduce the discrete version of the deviation of the total energy from the thermal equi- 
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librium for n > 0, 



£" = HK) (H(NZ) H{N%) - h{N%) (N% ~ N%)) 

KET 

+ ™( K ) ( H (Pk) H(P^) - &(P£)(Pg - P%)) 



KeT 



<j=K\L 



\ E E 

KGTcrESV 



and the discrete version of the energy dissipation ([7]): for n > 0, 



X" = J2 T ° min (N£ +1 ,N2 +1 ) \D (h (N n+1 ) - V n ] 



K.a 



aE£ int 
a=K\L 



■J2 E r a min (N-+\N^)[D(h(N^)-V n : 

KET aEC^t.K 

J2 T a wm{Pp 1 ,F£ +1 )[D(h(P»+ 1 )+V») Kt " 2 

aE£ int 
a=K\L 

E r CT min(P»+ 1 ,Pr 1 )[^(/ l (P"+ 1 )+^ 

KET aESevt.K 



K.a 



'K.a 



We present a test case for a geometry corresponding to a PN-junction in 2D picked in the paper of 
C. Chainais-Hillairet and F. Filbet [5]. The doping profile is piecewise constant, equal to +1 in the 
N-region and —1 in the P-region. 
The Dirichlet boundary conditions are 

77 = 0.1, P = 0.9, V 
77 = 0.9, P = 0.1, V 
Elsewhere, we put homogeneous Neumann boundary conditions. 

The pressure is nonlinear: r(s) — s 7 with 7 = 5/3, which corresponds to the isentropic model. 
We compute the numerical approximation of the thermal equilibrium and of the transient drift- 
diffusion system on a mesh made of 896 triangles, with time step At = 0.01. 

We then compare the large time behavior of approximate solutions obtained with the three following 
fluxes: 

— the upwind flux defined by (|2"31 (Upwind), 

— the Scharfetter-Gummel extended flux (l3ll) with the first choice (|27l) of drx.a, close to that of 
Jiingel and Pietra (SG-JP), 

— the Scharfetter-Gummel extended flux (f3~Tj) with the new definition (|2"9")l of dr^.a (SG-ext). 



h{N) - h{P) 



on {y = 1, < x < 0.25}, 



h(N) - h(P) 



on {y = 0}. 
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Fig. 3 Evolution of the relative energy £ n and its dissipation X n in log-scale for different schemes. 
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Fig. 4 The relative energy £ n and its dissipation I n in log-scale for different time steps. 



In Figure [3] we compare the discrete relative energy £ n and its dissipation I n obtained with the 
Upwind flux, the SG-JP flux and the SG-ext flux. With the third scheme, we observe that £ n 
and I n converge to zero when time goes to infinity, without a saturation phenomenon. This scheme 
is the only one of the three which preserves thermal equilibrium, so it appears that this property is 
crucial to have a good asymptotic behavior. 

In Figure 0] we compare the relative energy £ n and its dissipation X n obtained with the SG-ext 
flux for three different time steps At = 5.1CP 3 , 10~ 3 , 1CP 4 . It appears that the decay rate does not 
depend on the time step. 



5.2.2 The porous media equation 



We recall that the unique stationary solution u eq of the porous media equation (fTU|) is given by the 
Barenblatt-Pattle type formula (fTTj) , where Cj is such that u eq as the same mass as the initial data 
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uq. We define an approximation {Up-) KeT °^ u&q 



Tjeq 
K 



2l 



\xk\ 



1/(7-1) 



KeT, 



where C\ is such that the discrete mass of (U^) Ke7 - is equal to that of \Uk)kgT'> namer y 



E 



KeT KeT 

We introduce the discrete version of the relative entropy 



m(K)U^. We use a fixed point algorithm to compute this constant C\. 

12) 



KeT ^ 

and the discrete version of the entropy dissipation (|13l) 



Xk\ 



iU£-U%) 



E 

a=K\L 



T ff min ([/£,[/£) 



D h{U n ) + 



FT 
2 



RT,<r 



We consider the following two dimensional test case: r(s) 



D h(U n ) + 



2 




e-(x-2) 2 -(y+2) 2 

1 

' S-(x+2y-{y-2)* 






s , with initial condition 

{x - 2) 2 + (y + 2) 2 < 6, 
{x + 2) 2 + (y - 2) 2 < 6, 



and periodic boundary conditions. 

Then we compute the approximate solution on fl x (0,10) with Q = (—10,10) x (—10,10). We 
consider a uniform cartesian grid with 100 x 100 points and the time step is fixed to At = 5.10 -4 . 
In Figure we plot the evolution of the numerical solution u computed with the SG-ext flux 
at three different times t = 0, t = 0.4 and t = 4 and the approximation of the Barenblatt-Pattle 
solution. In Figure [5] we compare the relative entropy £ n and its dissipation T n computed with 
the scheme (|22| and different fluxes: the Upwind flux, the SG-JP flux and the SG-ext flux. We 
made the same findings as in the case of the drift-diffusion system for semiconductors: the third 
scheme is the only one of the three for which there is no saturation phenomenon, which confirms 
the importance of preserving the equilibrium to obtain a consistent asymptotic behavior of the 
approximate solution. Moreover it appears that the entropy decays exponentially fast, which has 
been proved in [7]. 

In Figure [7j we represent the discrete L 1 norm of U — U eq (obtained with the SG-ext flux) in log 
scale. According to the paper of J. A. Carrillo and G. Toscani, there exists a constant C > such 
that, in this case, 

\\u(t,x)-u eq (x)\\ L i {m) < Crap (-ft) , t >0. 



We observe that the experimental decay of u towards the steady state u eq is exponential, at a rate 
3 

better than -. 
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Fig. 5 Evolution of the density of the gas u and stationary solution u e 
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Fig. 6 Evolution of the relative entropy £ n and its dissipation I" in log-scale for different schemes. 



6 Conclusion 



In this article, we presented how to build a new finite volume scheme for nonlinear convection- 
diffusion equations. To this end, we have to adapt the Scharfetter-Gummel scheme, in such way 
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Fig. 7 Decay rate of \\U - U eq \\ L i. 



that ensures that a particular type of steady-state is preserved. Moreover, this new scheme is easier 
to implement than existing schemes preserving steady-state. 

In addition, we have shown that there is convergence of our scheme in the nondegenerate case. The 
proof of this convergence is essentially based on a discrete L 2 (0,T;if x ) estimate ([3"6|). A first step 
to then prove the convergence in the degenerate case would be to show this estimate without using 
the uniform lower bound of us- 

Finally, we have observed that this scheme appears to be more accurate than the upwind one, even 
in the degenerate case. Indeed, we have applied it to the drift-diffusion model for semiconductors 
as well as to the porous media equation. In these two specific cases, we clearly underlined the effi- 
ciency of our scheme in order to preserve long-time behavior of the solutions. At this point, it still 
remains to prove rigorously this asymptotic behavior, by showing a similar estimate to the one of 
the continuous framework (|5]) for discrete energy and discrete dissipation. 
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